In [1]:
"""
This first cell in the Jupyter Notebook is one long pseudo-comment
that shows my old solution set in Matlab syntax!
>> % a percent sign is how you write a comment in Matlab/Octave,
>> % so any line here that starts with a percent will be ignored by Matlab.
>> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> %%%%%%%%%%%%%%%%
>> % Question 1a
>> c=[3;4]
c =
3
4
>> x=[5;6]
x =
5
6
>> c' * x
ans =
39
>> %%%%%%%%%%%%%%%%
>> % Question 1b
>> % hard to do symbolic calculation in Matlab.
>> % Answer is 3*x_1 + 4*x_2
>> %%%%%%%%%%%%%%%%
>> % Question 1c
>> % c and x as above
>> c * x'
ans =
15 18
20 24
>> % notice that row 1 is 3*[5 6], and row 2 is 4*[5 6]
>> %%%%%%%%%%%%%%%%
>> % Question 1d
>> c=[3;4]
c =
3
4
>> B=[1,2;1,4]
B =
1 2
1 4
>> b=[6;12]
b =
6
12
>> c' * inv(B) * b
ans =
12
>> % just to see the middle step,
>> inv(B)
ans =
2 -1
-0.5 0.5
>> % or you could also do
>> B^-1
ans =
2 -1
-0.5 0.5
>>
>>
>> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> % Question 2: how many solutions to these systems?
>> % give one solution if one exists.
>> % if more than one solution exists, give at least two solutions.
>> %%%%%%%%%%%%%%%%
>> % Question 2a
>> A=[2,2;3,4]
A =
2 2
3 4
>> b=[5;6]
b =
5
6
>> inv(A)
ans =
2 -1
-1.5 1
>> % A is invertible (full rank), so there should be exactly one solution:
>> x = inv(A) * b
x =
4
-1.5
>> %%%%%%%%%%%%%%%%
>> % Question 2b
>> A=[1,2,3;4,5,6;7,8,9]
A =
1 2 3
4 5 6
7 8 9
>> b=[10;11;12]
b =
10
11
12
>> % is A invertible?
>> inv(A)
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.541976e-018.
ans =
-4.5036e+015 9.0072e+015 -4.5036e+015
9.0072e+015 -1.8014e+016 9.0072e+015
-4.5036e+015 9.0072e+015 -4.5036e+015
>> % Whoah, that doesn't look good. Let's try a different way:
>> rank(A)
ans =
2
>> % rank 2 for a 3-by-3 matrix means it's not invertible.
>> % So there should be either 0 solutions or infinite # of solutions.
>> % how to tell? Look at the "augmented" matrix:
>> [A, b]
ans =
1 2 3 10
4 5 6 11
7 8 9 12
>> % and look at the row-reduced form:
>> rref( [A, b] )
ans =
1 0 -1 -9.3333
0 1 2 9.6667
0 0 0 0
>> % Since that last row is 0 x_1 + 0 x_2 + 0 x_3 = 0,
>> % it's always satisfied, so the equations are consistent.
>> % Thus, there is an infinite # of solutions.
>> % How could we find them?
>> % or at least two of them?
>> % We could manipulate the RREF results above to get
>> % x1 = -9.333 + 1 x_3
>> % x2 = -9.6667 - 2 x_3
>>
>>
>> %%%%%%%%%%%%%%%%
>> % Question 2c
>> A=[0,2,3;4,5,6;7,8,9]
A =
0 2 3
4 5 6
7 8 9
>> b
b =
10
11
12
>> inv(A)
ans =
-1 2 -1
2 -7 4
-1 4.6667 -2.6667
>> % that looks nicely invertible.
>> x=inv(A) * b
x =
-5.3291e-015
-9
9.3333
>> % that first part of x is probably just 0.
>> %%%%%%%%%%%%%%%%
>> % Question 2d
>> A=[1 2 3; 4 5 6]
A =
1 2 3
4 5 6
>> b=[7;8]
b =
7
8
>> % A is clearly not invertible, since it's not square.
>> % More columns (variables) than rows (equations),
>> % so it should have an infinite # of solutions
>> % (unless it's the rare occurence of being inconsistent!)
>> % We could do RREF on the augmented matrix:
>> augmented = [A,b]
augmented =
1 2 3 7
4 5 6 8
>> rref(augmented)
ans =
1 0 -1 -6.3333
0 1 2 6.6667
>> % so then
>> % x_1 = -6.333 + 1 x_3
>> % x_2 = +6.6667 - 2 x_3
>> % Or, try this way:
>> A12 = A(:,1:2)
A12 =
1 2
4 5
>> A3 = A(:,3)
A3 =
3
6
>> % now pick an x3 and solve A12*[x1;x2] + A3*x3 = b
>> x3 = 0
x3 =
0
>> inv(A12)*(b-A3*x3)
ans =
-6.3333
6.6667
>> % and try a different x3
>> x3 = 1
x3 =
1
>> inv(A12)*(b-A3*x3)
ans =
-5.3333
4.6667
>>
>>
>> %%%%%%%%%%%%%%%%
>> % Question 2e
>> A=[1,2;3,4;5,6]
A =
1 2
3 4
5 6
>> b=[7;8;9]
b =
7
8
9
>> % the system is overdetermined (2 vars, 3 eqns) so unless we're really
>> % lucky, it will be inconsistent: no solutions.
>> rref([A,b])
ans =
1 0 -6
0 1 6.5
0 0 0
>> % This is saying
>> x=[-6; 6.5; 0]
x =
-6
6.5
0
>> A*x
??? Error using ==> mtimes
Inner matrix dimensions must agree.
>> x=[-6; 6.5]
x =
-6
6.5
>> A*x
ans =
7
8
9
>> % Wow, it isn't inconsistent! That's rare.
>> %%%%%%%%%%%%%%%%
>> % Question 2f
>> A=[3,5;21,35]
A =
3 5
21 35
>> b=[49;28]
b =
49
28
>> x=inv(A)*b
Warning: Matrix is singular to working precision.
x =
Inf
Inf
>> % well that's bad news. Let's check the rank:
>> rank(A)
ans =
1
>> % Ah, the second rows is 7 times the first row
>> % (or, the 2nd column is 5/3 times the first column!)
>> % but 28 is not 7 times 49.
>>
>>
>>
>>
>> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% switched to Octave at this point.
>> % Question 4: Graph the set of solutions to
> [x,y]=meshgrid(linspace(0,10,20),linspace(0,10,20));
> z = 1*x + 2*y;
> contour(x,y,z,[10 10])
> % the [10 10] is a little strange;
> % We want just the z=10 contour, but if we say
> % contour(x,y,z,10)
> % it will interpret it as "give me 10 contour lines, please"
> % so we have to pretend to give it a list of the contour values we want.
> % (a list being at least two elements)
>
>> %%%%%%%%%%%%%%%%
> % Question 4b: [3 4][x;y] <= 12
> z=3*x+4*y;
> feas = z <= 12;
> contour(x,y,feas,[0.5 0.5])
> % ew, that's ugly. Too jagged.
> % Could do it this way:
> contour(x,y,z,0:12)
> % to show all the contours 0,1,2,...,12
> % or, do the "feas" thing but with better resolution
> [x,y]=meshgrid(linspace(0,10,50),linspace(0,10,50));
> z=3*x+4*y;
> feas = z <= 12;
> contour(x,y,feas,[0.5 0.5])
>> %%%%%%%%%%%%%%%%
> % Question 4c: [3 4;6 2][x y] = [12;12]
> z1=3*x+4*y;
> z2=6*x+2*y;
> contour(x,y,z1,[12 12])
> hold on
> contour(x,y,z2,[12 12])
> A=[3 4;6 2];
> xvec=[4/3 , 2]
xvec =
1.3333 2.0000
> A*xvec
error: operator *: nonconformant arguments (op1 is 2x2, op2 is 1x2)
error: evaluating binary operator `*' near line 25, column 2
> A*xvec'
ans =
12
12
> plot(4/3, 2.0, 'o')
"""
0
Out[1]:
In [2]:
# a number-sign (a.k.a "hashtag" these days) is how you write a comment in your Python code:
# anything after that (whether it starts the line or comes after useful code on the line)
# will be ignored.
# first, we'll import the usual libraries that we need:
import numpy as np
import matplotlib.pyplot as plt
# and do the magic to make plots appear in the notebook:
%matplotlib notebook
In [3]:
# There are various ways to define things that are like vectors or matrices,
# and not all of them are very helpful.
# If you're used to Matlab you'd do this for a row vector:
c = [3,4]
c
Out[3]:
In [4]:
# and you get back what you typed, but it's actually a python "list":
type(c)
Out[4]:
In [5]:
# lists are useful too, but they aren't the best to do linear algebra with.
# instead we need a NumPy "array", which is also called an ndarray (stands for n-dimensional array)
# Note that there's also a NumPy "matrix" class but it's less useful.
# And, there's a plain-python "array" class but again it's less useful.
c=np.array([3,4])
c
Out[5]:
In [6]:
type(c)
Out[6]:
In [7]:
# how big is it? rows? columns?
# the "shape" attribute tells us (similar to matlab's "size" function,
# but here we don't need parentheses like c.shape() because .shape is not a function )
c.shape
Out[7]:
In [8]:
# saying (2,) is python's way of saying it's 1-dimensional, so it really
# only has a length rather than length-by-width. It's not specifically
# a row and not specifically a column, it's just a vector.
# The comma with nothing after it is python's way of saying it's
# like an ordered pair or ordered triplet, but with only 1 thing: an ordered single?
# The generic name for ordered pairs, ordered triplets,
# ordered quadruplets/quintuplets/sextuplet, etc. is a "tuple".
# see https://docs.scipy.org/doc/numpy-dev/user/quickstart.html
# If you want to force it to be 2-by-1 or 1-by-2 you can do it this way:
c=np.array([ [3,4] ]) # a row vector, 1-by-2; note the extra set of [] brackets
c
Out[8]:
In [9]:
c.shape
Out[9]:
In [10]:
# or this way:
c=np.array([ [3], [4] ]) # a column vector, 2-by-1
c
Out[10]:
In [11]:
c.shape
Out[11]:
In [12]:
# but putting the brackets around each component can get tiring.
# Instead, you can create it as a row vector and immediately
# transpose it:
c=np.array([ [3,4] ]).T # a row vector, 1-by-2, then the .T transposes it:
c.shape
Out[12]:
In [13]:
# I think the usual practice in Python is to treat a vector
# as a 1-dimensional thing like c=np.array([3,4]) with shape (2,)
# but since we're just learning, and we know how linear algebra
# works with row vectors distinct from column vectors,
# I'm going to use things that are clearly a row or column vector.
In [14]:
# Question 1a
c=np.array([[3,4]]).T # ends up as a column vector
x=np.array([[5,6]]).T # ends up as a column vector
c.T @ x # the @ is how you do matrix-multiplication in Python/numpy.
Out[14]:
In [15]:
# Whoah, that's interesting--it returns a 1-by-1 array rather than just calling it a scalar.
# If you think of defining the dot-product a dot b as the sum of (a_i)*(b_i), that returns a scalar
# (not a 1-by-1 matrix) and in fact the dot product is sometimes called the scalar product.
# But we usually think of dot product as the same thing as (vectortranspose)*othervector.
# The math community hasn't even really decided which one is preferable.
# We'll just move on.
In [16]:
#Question 1b
# It's hard to do symbolic calculation in Python, though you could use the
# sympy package (stands for symbolic-math)
# note that sympy is different than simpy, which does simulation of things like queues.
# Answer is 3*x_1 + 4*x_2
In [17]:
# Question 1c
# c and x as above
c @ x.T
#ans =
# 15 18
# 20 24
#notice that row 1 is 3*[5 6], and row 2 is 4*[5 6]
Out[17]:
In [18]:
#Question 1d
c=np.array([[3,4]]).T # ends up as a column vector
B=np.array([[1,2],[1,4]])
b=np.array([[6,12]]).T # ends up as a column vector
c.T @ np.linalg.inv(B) @ b
Out[18]:
In [19]:
# just to see the middle step there:
np.linalg.inv(B)
Out[19]:
In [20]:
# Note that if you try to "raise B to the -1 power" this way:
B**-1.0
# you get 1/(each element of B individually), not what we wanted!
Out[20]:
In [21]:
###############################
# Question 2: how many solutions to these systems, A*x=b?
# If a solution exists and is unique, give it.
# If more than one solution exists, give at least two solutions.
###########
# Question 2a
A=np.array([[2,2],[3,4]])
b=np.array([[5,6]]).T # ends up as a column vector
# it's 2 equations in 2 variables, so it has some hope of a unique solution.
# We could check if A is invertible.
np.linalg.inv(A)
Out[21]:
In [22]:
# A is invertible (full rank), so there should be exactly one solution
# We could find it via
np.linalg.inv(A) @ b
Out[22]:
In [23]:
# but taking the matrix inverse is slow compared to solving for just one "x" vector value,
# using Gaussian Elimination or something like that (Gaussian Elimination itself isn't
# all that great in many cases; the .solve function uses something fancier):
np.linalg.solve(A,b)
Out[23]:
In [24]:
#Question 2b
A=np.array([[1,2,3],[4,5,6],[7,8,9]])
b=np.array([[10,11,12]]).T # ends up as a column vector
# it's 3 equations in 3 variables, so it has some hope of a unique solution.
# We could check if A is invertible:
np.linalg.inv(A)
Out[24]:
In [25]:
# Hmm, it didn't complain but we got a strange matrix with huge values,
# which could have been caused by division-by-near-zero.
# One way to think computing the matrix inverse involves dividing by its determinant;
# if that was zero or near-zero we could have trouble. Let's check:
np.linalg.det(A)
Out[25]:
In [26]:
# that's essentially 0, so the A matrix is not invertible. We could also ask its rank.
# If the rank is less than 3 for this 3-by-3 matrix, the matrix isn't invertible.
np.linalg.matrix_rank(A)
Out[26]:
In [27]:
# We could still see what we get if we try to solve:
x=np.linalg.solve(A,b) # solve A*x=b to get x
x
Out[27]:
In [28]:
# it gave us a solution, and didn't complain. Does that mean
# everything is going well? We could at least check the solution:
A @ x
# compare the result to b=[10,11,12]
Out[28]:
In [29]:
# So at least one solution exists, and we know the matrix is non-invertible,
# so there must be other solutions.
# We could get them by forming the augmented matrix [A,b] and doing RREF
# on it, then figuring out the general form of the solution.
# It turns out Python doesn't have a built-in RREF in NumPy (it does in SymPy)
# The RREF turns out to be
# 1 0 -1 -9.3333
# 0 1 2 9.6667
# 0 0 0 0
#Since that last row is 0 x_1 + 0 x_2 + 0 x_3 = 0,
#it's always satisfied, so the equations are consistent.
# Thus, there is an infinite # of solutions.
# How could we find them?
# or at least two of them?
# We could manipulate the RREF results above to get
# x1 = -9.333 + 1 x_3
# x2 = -9.6667 - 2 x_3
In [30]:
# Question 2c
A=np.array([[0,2,3],[4,5,6],[7,8,9]])
b=np.array([[10,11,12]]).T # ends up as a column vector
# it's 3 equations in 3 variables, so it has some hope of a unique solution.
# We could check if A is invertible:
np.linalg.inv(A)
Out[30]:
In [31]:
# Looks nicely invertible. Now let's solve by doing inv(A)*b:
np.linalg.inv(A) @ b
Out[31]:
In [32]:
# that first component is probably really just 0.
In [33]:
# Question 2d
A=np.array([[1,2,3],[4,5,6]])
b=np.array([[7,8]]).T # ends up as a column vector
# A is clearly not invertible, since it's not square.
# It has more columns (variables) than rows (equations),
# so it should have an infinite # of solutions
# (unless it's the rare occurence of being inconsistent!)
# We could do RREF on the augmented matrix and get
# 1 0 -1 -6.3333
# 0 1 2 6.6667
# so then
# x_1 = -6.333 + 1 x_3
# x_2 = +6.6667 - 2 x_3
# But, here's a completely different way that is important later in Math 560.
# If we have freedom in 1 of our 3 variables (could be any one),
# let's just pick x3. We could set it to anything we want.
# Zero is a very nice thing to set it to.
# Then we'll rewrite the system with our choice of x_3 on the right-hand side,
# no longer acting as a variable, and get a 2-by-2 system.
# First we'll take the 3rd column of A by itself:
Acol = A[:,2] # this is saying: all rows (the : means "all"), and column#2
Acol
# NOTE!!!!!!!!!!!!!!!!!
# Python uses column (and row) numbers that start at 0, not 1 !!!!!!!!!!!
# so if we want what linear algebra people call the 3rd column, we say column 2.
Out[33]:
In [34]:
# Hm, it ended up as a 1-dimensional vector (neither specifically row or column),
# rather than as a 2-by-1 (column) matrix.
#Then we look at the rest of A, which is square.
# NOTE!!!!!!!!!!!!!!!!!
# The indexing here gets even stranger: saying 0:5, for example, produces 0,1,2,3,4 BUT NOT 5!!!!
# So if we want columns 0 and 1, we say 0:2
Asquare = A[:,0:2]
Asquare
Out[34]:
In [35]:
# Now we solve A*x = b by breaking it up:
# Asquare * [x1;x2] + Acol * x3 = b
# then bring the x3 to the other side:
# Asquare * [x1;x2] = b - Acol * x3
# but substitute in whatever value of x3 we want:
x3 = 0
np.linalg.solve(Asquare, b - Acol*x3)
Out[35]:
In [36]:
# and do it again for any other value of x3 we want:
x3 = 1
np.linalg.solve(Asquare, b - Acol*x3)
Out[36]:
In [37]:
# If you did this via matrix-inverse and used the same Asquare I did, you used
np.linalg.inv(Asquare)
Out[37]:
In [38]:
#Question 2e
#A=[1,2;3,4;5,6]
A=np.array([[1,2],[3,4],[5,6]])
b=np.array([[7,8,9]]).T # ends up as a column vector
# the system is overdetermined (2 vars, 3 eqns) so unless we're really lucky,
# it will be inconsistent: no solutions.
# We could form the augmented matrix [A,b] and finds its RREF, and get
# 1 0 -6
# 0 1 6.5
# 0 0 0
# notice that last line says 0*x1+0*x2 = 0, which is consistent!
# The other two rows say 1*x1 + 0*x2 = -6 or x1=-6
# and 0*x1+1*x2 = 6.5 or x2=6.5, so that's our solution: [-6;6.5]
# If we hadn't done the RREF, we could just pick two rows from A
# instead of all 3, and solve and hope for the best:
whichrows = range(0,2)
Asquare = A[whichrows,:]
bsmall = b[whichrows,:]
np.linalg.solve(Asquare,bsmall)
Out[38]:
In [39]:
# Because the system is consistent, we could try the last 2 rows instead:
whichrows = range(1,3)
Asquare = A[whichrows,:]
bsmall = b[whichrows,:]
np.linalg.solve(Asquare,bsmall)
Out[39]:
In [40]:
# or the first and last row:
whichrows = [0,2]
Asquare = A[whichrows,:]
bsmall = b[whichrows,:]
np.linalg.solve(Asquare,bsmall)
Out[40]:
In [41]:
# a more general solution for an overdetermined system,
# in case it's not consistent, is to ask for a solution
# that makes A*x as close as possible to the right-hand-side b,
# in a least-squares sense:
x, resid, rank, singularvalues = np.linalg.lstsq(A,b)
x
# This is actually equivalent to doing linear regression, solving
# (first column of A) * slope1 + (second column of A)*slope2 = yvalues (or bvalues, same thing)
Out[41]:
In [42]:
# Question 2f
A=np.array([[3,5],[21,35]])
b=np.array([[49,28]]).T # ends up as a column vector
np.linalg.solve(A,b)
In [43]:
# It says the matrix is singular.
# Ah, the second row is 7 times the first row
# (or, the 2nd column is 5/3 times the first column!)
# but on the RHS, 28 is not 7 times 49.
# so the system is inconsistent; there are no solutions.
In [44]:
# Question 4a: Graph the set of solutions to
# [1 2][x; y] = 10
# You could do this by going to Desmos.com and typing
# 1*x + 2*y = 10
In [45]:
# Here's a crazy/good/crazy-good way to do it: think of
# 1*x + 2*y = 10
# as graphing the value-10 contour of the function f(x,y)=1*x+2*y
import matplotlib
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
matplotlib.rcParams['xtick.direction'] = 'out'
matplotlib.rcParams['ytick.direction'] = 'out'
delta = 0.5
x = np.arange(-1.0, 10.0, delta)
y = np.arange(-1.0, 10.0, delta)
X, Y = np.meshgrid(x, y)
Z = 1*X + 2*Y
plt.figure()
levels = [10]
CS = plt.contour(X, Y, Z,levels)
plt.clabel(CS, inline=1, fontsize=10)
plt.title('Simplest default with labels')
Out[45]:
In [46]:
# Question 4b: Graph the set of solutions to
# [3 4][x; y] <= 12
# You could do this by going to Desmos.com and typing
# 3*x + 4*y <= 12
# and it adds a color to the side that is feasible.
In [47]:
# Question 4c: Graph the set of solutions to
# [3 4;6 2][x; y] = [12;12]
# You could do this by going to Desmos.com and typing
# 3*x + 4*y = 12
# 6*x + 2*y = 12
# It will show 2 lines and their intersection.
# If you click on the intersection, it will show its coordinates.
# Technically, if we're just graphing the "set of solutions" to the system,
# the only thing to graph is that 1 point that is the intersection.
# But it's nice to show what goes into determining that.
# We can also find that point by using linear algebra:
A=np.array([[3,4],[6,2]])
b=np.array([[12,12]]).T # ends up as a column vector
np.linalg.solve(A,b)
Out[47]:
In [ ]: